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Brownian dynamics algorithms integrate numerically Langevin equations and allow to probe long 
time scales in simulations. A common requirement for such algorithms is that interactions in the sys- 
tem should vary little during an integration time step: therefore, computational efficiency worsens 
as the interactions become steeper. In the extreme case of hard-body interactions, standard nu- 
merical integrators become ill defined. Several approximate schemes have been invented to handle 
such cases with little emphasis on testing the correctness of the integration scheme. Starting from 
the two-body Smoluchowsky equation, we discuss a general method for the overdamped Brownian 
dynamics of hard-spheres, recently developed by one of us. We test the accuracy of the algorithm 
with the exact solution of the Smoluchowsky equation in the case of two body collisions and in the 
low-density limit. 

PACS numbers: 05.40.Jc, 05.10.Gg, 61.20.Ja 



I. INTRODUCTION 

The simulation of interacting Brownian particles, 
called Brownian Dynamics (BD) simulation [ij, has be- 
come an important tool in condensed matter, colloidal 
and biological physics. In particular at high densities, 
the excluded volume has a major effect on the dynamics 
of the Brownian particles. Consequently, many model 
colloids are well described by effective pair interactions 
that are hard-sphere (HS) like, or have a HS core with ad- 
ditional potential tails. This places the hard-sphere sys- 
tem (with or without random 'Brownian' forces) among 
the most important reference systems for the theories in 
the field 0, Q- However, such theories are in general 
only approximate, requiring extensive testing through 
experiment and computer simulation. Yet, the simu- 
lation of hard spheres with Brownian dynamics is less 
than straightforward, because of the singular nature of 
the interaction potential: most methods dealing with the 
numerical integration of stochastic differential equations 
(SDEs) require a certain degree of smoothness in all the 
interactions. 

Up to now, several hard-core BD algorithms have been 
proposed and apphed in the past with 
a varying degree of justification. Our aim here is to re- 
view these approaches and to compare them with a novel 
approach developed by one of us [ijl, called De Michele's 
algorithm in the following. We investigate the conver- 
gence to the true solutions of the singular SDE by study- 
ing special situations where the theoretical behavior is 
still under control, in order to validate the correctness of 
the approach. This is the first step in a programme to ex- 
tend De Michele's algorithm to more complicated cases, 
such as taking into account inertial effects (assumed to be 
negligible in proper BD), or more complicated singular 
interactions. 

Of course the hard-sphere potential (/31/(r) = cx) if two 



particles overlap, zero otherwise, where /3 = l/(fcBT) is 
the inverse temperature), is a purely theoretical concept. 
One may view steeply repulsive 'soft-sphere' potentials, 
e.g., V{r) (X r~" with large n, as more convenient model 
systems for colloids |5lll2j and even try to infer HS (n 
oo) behavior by mapping the n-dependent results, taking 
into account structural information about the system 

. Such soft-sphere methods are of course hampered by 
having to use smaller and smaller integration time-steps 
for increasing n and careful extrapolation to n — > cx). 
Therefore, a proper HS-BD algorithm is worthwhile. 

We will not discuss hydrodynamic interactions (HI), 
i.e., the solvent-induced interactions that are present 
in typical experimental realizations of Brownian sys- 
tems. Taking them into account properly ensures that 
hard s phe res cannot overlap, due to divergent lubrication 
forces |lj] . To deal with HI, several computational meth- 
ods have been developed, for example Stokesian Dynam- 
ics ^3 1 Lattice Boltzmann simulations [l^ . Il^ lisl , 
Dissipative Particle Dynamics [20I l2l| , or Fluid Particle 
methods [23. Usually they either deal with softened in- 
teractions again and/or are rather time-consuming: de- 
pending on the method, a huge number of degrees of 
freedom needs to be considered, or the intricate nature 
of the non-pairwise-additive long-range HI forces the use 
of elaborate schemes. While the theory of HI is rather 
well understood at low particle densities, much less is 
known at high densities, and theories often proceed by 
claiming them irrelevant. A non-HI simulation therefore 
still has its place in testing such theories, and in circum- 
venting the huge effort of the HI methods, should the 
claim be true. 

Testing the goodness of the approximations inherent 
to all existing HS-BD algorithms has up to now received 
little attention. One needs to test properties that are in- 
herent both to the Brownian dynamics of the system and 
to the hard-core collisions. Testing for diffusive behavior 
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is obviously not enough, in as far as at high densities and 
long times it is related to the chaotic nature of the many- 
body problem and not to the details of the implemented 
equation of motion or its correctness. In some previous 
work, static quantities like the radial distribution func- 
tion g{r) have been checked 0, but such comparisons 
do not test whether an HS-BD algorithm properly dis- 
cretizes the Langcvin SDE, but rather its crgodicity. 

We will therefore discuss tests of BD-HS algorithms 
that involve probing the hard-sphere interactions but still 
allow for a comparison with exact results known for the 
Brownian system, i.e. where the Langevin equation can 
be solved exactly in terms of a time-dependent proba- 
bility distribution function (PDF). In such a case, an 
empirical PDF can be generated from sufficiently many 
runs of the algorithm in question and then be compared 
to the exact solution. 

The methods we are going to discuss in the present 
paper, are crossbreeds between standard Brownian meth- 
ods (ignoring the singular potential), and the standard 
simulation method for the ordinary (non-stochastic) dy- 
namics of hard spheres, called event-driven (ED) simu- 
lations. Predecessors date back to Monte Carlo (MC) 
inspired schemes 0, S 13 1 culminating in the algo- 
rithm by Strating and algorithms making use of so- 
called overlap potentials Strating's algorithm is but 
a step away from De Michele's algorithm [llj. A slightly 
different principle has also recently been implemented by 
Tao et al. |lCll |. albeit applied to a more complex system 
of hard rods. We will discuss both these event-driven 
Brownian dynamics (ED-BD) methods, with an empha- 
sis on De Michele's algorithm, which we will show to be 
exact to first order in the integration-step size At. 

This paper is organized as follows: the next sec- 
tion discusses the theoretical background, introduces 
De Michele's algorithm and extends it to the case of shear 
and constant forces. Section IlIII presents numerical tests 
for two-particle cases where a comparison with the exact 
PDFs is possible. In section ITVI we compare the diffusive 
many-body behavior with exact theoretical predictions. 
Section El summarizes and concludes. 

II. HARD-SPHERE BROWNIAN DYNAMICS 
ALGORITHMS 

We start from the iV-particle (1 < i < N) Langevin 
equation 

mA^n+f\+fi+fr, (1) 

where the subscripts i label the particles and will be 
dropped where they are clear from the context. In 
Eq. J^l, incorporates the effect of the HS interac- 
tion and /°'^' is a sufficiently smooth external force. /' 
and Z*^ are the random and dissipative forces resulting 
in Brownian dynamics, therefore the above equation is a 
second-order SDE in the particle positions. The random 



force is characterized by its zero mean and the correlation 

{T{t)®f{t'))^2k^TR5{t-t'), (2) 

and the fluctuation-dissipation theorem (FDT) then re- 
quires the dissipation to be j = —R{v — u), where u is a 
local flow velocity field; u is nonzero in the case of applied 
shear, u = Fr, where F is a shear-rate tensor. R is in 
general a complicated matrix depending on the full con- 
figuration of the system at any given time, representing 
hydrodynamic interactions (HI). Here we are concerned 
with the simpler case where i? = ^1 is constant, with a 
real number ^ > characterizing the noise. In this paper 
we furthermore deal with the limit m/^ — > 0, the case of 
strong dissipation, where Eq. reduces to 

s,v^n+fi+iu+fT\ (3) 

that is a first-order SDE in the position of the particles. 

In what follows, we will, according to customary pro- 
cedure |23| , assume that one can fix a small enough time 
interval At over which the particle configuration and all 
smooth forces vary slowly, allowing them to be treated as 
constant. If this were possible for as well, one could 
use conventional SDE integrators [2 51 . a topic in its own 
right (see e.g. Refs. i ^ ^ iMlS, El 111 111 111 ) . 
For hard spheres, the forces are not Lipschitz continuous 
and standard numerical integrators are ill-defined |35l |. 
Integrating Eq. ^ over the interval At, one obtains 

Af(At)=5At + M(At), (4) 

where g = u + ^~^{f^ + /°'^') contains the systematic 
and interaction terms, while /i is a Wiener process with 

(/i(t) (g) fl{t')) = 2Dmm{t,t'), where D = k^TR^^^ = 
(/cbT/^)! is the matrix of diffusion coefficients. Note that 
we still allow these bare diffusion coefficients to depend 
on the particle index, as it is the case in multi-component 
(polydisperse) systems. 

We discuss the case u = f'^^^ = first. One approach 
is to simply set = in generating random displace- 
ments Ar with the statistics given by fl 0,0,0,0jI1jI3 ■ A 
second step then is needed to ensure that the random dis- 
placements are compatible with the presence of ^ 0; 
for hard spheres this is the criterion that no two par- 
ticles overlap. One can either discard all displacements 
that violate this condition - this is a variant of standard 
Metropolis MC (called simply MC in the following). It 
has been shown d that in the limit At one indeed 
obtains Brownian dynamics, i.e., a faithful realization of 
Eq. |(2J). We will discuss this point again further on. 

More elaborate overlap removal can be applied: some 
have put particles back at contact (or around contact 
on average) along the line of their relative displacement 
[1H0,|1|, although this has unwanted effects on the pair 
distribution function, hence on observablcs like pressure. 
Strating |1| proposed to remove overlaps by performing 
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elastic collisions: from the time step At and the displace- 
ment Ar , assign a fictive velocity v = Af/ At to the par- 
ticles, and move them by displacements Ar as if they 
underwent ballistic flight (including eventual elastic col- 
lisions) with this velocity v in the time interval At. 

Whatever the removal procedure, one has to be aware 
that the removal of one overlap can create another, called 
secondary (tertiary, and so forth) overlaps. This is espe- 
cially likely to happen in dense systems. As Strating 
pointed out, both the removal of all these 'higher-order' 
overlaps as well as the order in which overlaps are re- 
moved are crucial for the algorithm to work. According 
to Ref. 0, one needs to remove them in the order in 
which they would have occurred in a Newtonian dynam- 
ics simulation. 

At this point it is convenient to turn around the dis- 
cussion and start from an algorithm that treats Eq. (jS)) 
in the opposite limit: If ^ were zero, one would recover 
the standard Newtonian dynamics of hard spheres. The 
best numerical schemes in this case are event-driven (ED) 
simulations |36l l37j : assuming the collisions are binary 
in all non-degenerate cases and of infinitesimal duration, 
one advances the system from one such collision to the 
next, solving the ballistic free flight in between exactly. 
The random force in Eq. Q prevents a nai've application 
of this approach, but at least every At we can reintro- 
duce the randomness: if we interpret the velocities of 
the ED simulation as the fictive velocities of the above 
discussion, randomly drawn every At, the ED scheme is 
simply a device to prevent all unphysical overlaps in the 
first place, using a set of vectors v as its book-keeping 
device. 

Note that in terms of efficiency, there is little difference 
to Strating's method, since the main effort in ED simula- 
tion goes into the calculation and the sorting of collision 
times, something that is also needed when one is to re- 
move overlaps in their 'correct' order. Note also that 
this algorithm bears a resemblance to some ED granular 
matter simulations 39, 40,0,^2; where however a 
modified kinetic equation instead of Eq. ^ needs to be 
solved. 

The question then arises how to handle 'Brownian col- 
lisions' in such a combined ED-BD scheme. We will argue 
below that the elastic collision rule of Newtonian dynam- 
ics is a reasonable choice, inspired by its exactness in one 
dimension. This is De Michele's algorithm, first used 
in [Tlj . In the limit of small At and for vanishing u it 
is essentially Strating's algorithm, with the overlap prob- 
lem cured. (For large At, the two algorithms differ in the 
treatment of one particle 'tunneling' across another.) We 
we also discuss another choice, made by Tao et al. 
in the more complicated context of hard rods: instead of 
ballistic collisions, where the pre- and post-collision ve- 
locities are perfectly correlated, one can decorrelate them 
by assigning the post-collision velocity randomly (with 
proper restrictions to avoid overlaps). 

To investigate the role of ED collisions in the BD al- 
gorithm, we study the two-body BD-HS problem with- 



out external forces in more detail. This can be solved 
analytically [i^ . |43 | , by transforming from the particle 
coordinates r^, i = 1,2, to the relative, r = ri — r2, 
and center-of-motion coordinates, R = E!~^(^iri -l-^2'^2), 
where 5 = -I- ^2- The latter then separates from the 
problem, giving free center-of-motion diffusion. In the 
relative coordinates, the hard-core interactions can be 
included as a boundary problem into the corresponding 
diffusion equation, 

dtp{r, t) = DV^p{f, t), r > cr , (5a) 

r-Vp{f,t) = 0, r = <7, (5b) 

where p(r, t) is the PDF for finding a relative distance f 
at time t, D = Di + D2 is the relative diffusion coefficient 
and a = (cti -I- (T2)/2 with the particle diameters at. 

One could aim to include the exact solution of this 
case in an HS-BD algorithm, by drawing particle dis- 
placements according to this p{f, t) whenever two parti- 
cles are suflaciently near. However, the procedure would 
be rather involved, since the analytical solution 0,^3 
only available as an infinite sum in the Laplace domain. 
It would also be only approximate in the presence of more 
than two particles, hence valid only for small time steps, 
and again introduce secondary overlap problems. 

It is therefore easier to consider the limit of small At; 
then, the boundary at r = cr can be considered flat, and 
the analytical PDF is easily obtained by the method of 
images [43: denoting by Go(r, to + At|ro, to) the Green's 
function of the unbounded diffusion equation (for a par- 
ticle starting at tq at time to), we have 

p{f,to + At\ro,to) 

= Go(r, to + At|fo, to) + Go{r, to + At\r*,to) (6) 

for r > a, and zero otherwise. Here, ^o* the mirror 
image of the initial coordinate 7% with respect to the 
boundary. An algorithm which implements transitions 
of a particle from ro to f in the time interval At by 
performing a ballistic free flight with a velocity v cho- 
sen from some probability distribution f{v) that satisfies 
f{v = Af/At)dv = Go(Ar,to + At\0,to)dr, will satisfy 
the correct diffusion equation for such no-flux boundary 
conditions. This holds as long as the correct solution can 
be obtained by the image method, and if Gq describes 
a homogeneous process, i.e., depends on position only 
through the difference Ar = r — fo, and not on r and ro 
individually. 

The fictive-velocity distribution correspond- 
ing to free Brownian diffusion, Go(Ar, At) = 
exp[— (Ar')^/(4DAt)], is of course 
just the Maxwell-Boltzmann distribution, f{v) — 
(/3m/27r)~''/^exp[-/3(m/2)'(;^]. This connection is 
essentially the way in which the FDT enters the ED-BD 
algorithms. The parameter f]mi sets the bare diffusion 
coefficient for each particle i, given as 
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In particular, polydisperse systems, where different HS 
species liave different Di, can be simulated by incorpo- 
rating elastic collisions according to the (fictivc) masses 
rui. 

De Michele's algorithm for hard spheres then 
reads: (i) every t„ = nAt {n integer) extract veloci- 
ties Vi according to a Maxwellian distribution with fictive 
masses obeying Eq. (ii) evolve the system between 
tn and tn + Ai according to the laws of ballistic motion 
(performing ED molecular dynamics). 

As outlined above, one can anticipate that this algo- 
rithm converges to the correct solution in the limit At 
0, where the problem reduces to the one-dimensional 
two-body case, which is treated exactly. A first esti- 
mate for its regime of applicability is At <C (T^/(4£') 
- in order to treat particle boundaries as flat walls - and 
At ^ d^yg/(4D), where davg ~ p~^l^ is a typical inter- 
particle separation - so that 'almost all' time steps deal 
with binary interactions only. 

The image-method solution can be extended in a 
straightforward manner to include a linear shear field, 
u = Fr, and external forces that change slowly in space 
(such that they can be approximated as constant over 
the typical Ar resulting in a time step). Let us dis- 
cuss the case of linear shear flow in more detail. In this 
case one needs to deal with two effects on the Brow- 
nian motion: there will be a systematic drift u, but 
also the noise term for Ar will be modified because 
u depends on r. If the shear flow acts along the x- 
direction, u = jyes, we have 0: (Ay) = (Az) = 0, 
(Ax) = yjAt, and (Ay^) = (Az^) = 2D At, but (Ax^) = 
(Ax)^ + 2DAt{l - (1/3)(7A<)2). The random displace- 
ments appearing in Eq. Q in the presence of linear shear 
are also cross-correlated, {Ax Ay) = {DAt){'^At). 

This leads to the following extension of Do Michele's 
algorithm: (i) extract random velocities v = 
(Ar — (Ar)) / At from a multivariate Gaussian distribu- 
tion according to the above averages, (ii) add to the fic- 
tive velocities the systematic drift (Ar^ /At induced by 
the shear flow and/or gravity. In previous applications 
of HS-BD algorithms to sheared systems, only the sys- 
tematic drift has been taken into account (cf. Ref. |3|). 
The effect of neglecting these high-shear-rate corrections 
to the Brownian displacements has not been studied so 
far and remains to be clarified. For constant external 
forces (such as gravity), only a systematic drift needs to 
be taken into account. 



III. TWO-PARTICLE TESTS 

As in two (and higher) dimensions, no exact solution in 
terms of a mirror image exists for the two-body HS prob- 
lem, De Michele's algorithm introduces an approximation 
which is worthwhile testing. The exact distribution func- 
tion in two dimensions is. after a Laplace transform from 



time t to frequency s = D q |4 
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p{r,cb,s\ro,cl>o) = - ^ e™(*-*'')/(2^)if™(<jr>) 

, (8) 



/„,(gr<)-A™(gr<)|fM 

where we have transformed to polar coordinates (r, (jj) . 
Here, r> is the greater of r, rg, and r< the lesser. 
/m(z) and Km{z) arc the modified Bessel functions 
of the first and second kind of order m. Note 
that here and in the following, D denotes the rela- 
tive diffusion coefficient between the two spheres. On 
the other hand, the PDF of De Michele's algorithm, 
Eq. ©, is in two dimensions Laplace transformed us- 
ing LT[l/(47ri:)<)exp[-aV(4Dt)] = 1/ {2T:)KQ{aq). The 
latter may be expressed usin g th e addition theorem for 
the modified Bessel functions |43 as 
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, (9) 



where (rQ,(/)5) are the angular coordinates of the mirror 
image. Since they depend in a complicated fashion on 
(r, (j)) , it is hard to assess the error made in this approxi- 
mation analytically. At least in the limit q s- oo, one rec- 
ognizes from I^{qr*^)/Km{qr>) r^{qr%) / K'^{qr>), 

that J3ed('", 't'j s) approaches the true p{r, (f>, s) in the limit 
ro a. 

We can however proceed by looking at a particular 
average that can be calculated exactly. Let us con- 
sider the average displacement magnitude |(Af)(iD| = 
i^i{Axi)^j^y^^ where Xi label the Cartesian compo- 
nents and {Axi)dD ~ r'^^^ drd^dAxip{r, s|ro) is the 
d-dimensional average. With the above form of the 
two-dimensional exact PDF, Eq. ©, this average eval- 
uates to a relatively simple expression (see Appendix 
In three dimensions, the calculation proceeds along 
the same lines, essentially replacing the Bessel functions 
Im{z) and Km{z) by their spherical counterparts, im{z) 
and km{z), and involving spherical harmonics in place of 
the angular exponentials. We get 
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(11) 

(12) 



always using ro > cr. The one-dimensional result can also 
be obtained by recurring to the corresponding Laplace- 
transformed PDF or obtained in the time domain 
by direct integration (making use of the image-method 
solution). 
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FIG. 1: Average displacement size, |{Ar)2D|, for the two- 
dimensional Brownian motion of two equal hard spheres 
{a = 1) in their relative-coordinate frame; ro is the initial 
center-to-center distance. The relative diffusion coefficient is 
D = Af/2; the step size At is kept fixed within each panel. 
All values are scaled by 1/ \/DAt. Symbols are one-timestep 
simulation results, averaged over 50000 runs for each point: 
De Michele's algorithm (circles), the Tao-BD algorithm (plus 
symbols), a modified Tao-BD algorithm (crosses), Strating's 
algorithm (diamonds), and a MC algorithm (triangles). For 
At = 0.01, diamonds and circles, as well as crosses and plus 
symbols overlap. Solid lines are the analytical solutions. The 
dotted line corresponds to exp[— (ro — u)^/ (4Dt)]/\/ irDt. 



Note that the two- and three-dimensional expressions 
for large q have Eq. Hl()|l as their leading asymptote. The 
same holds for the PDF implemented by De Michele's 
algorithm, again demonstrating that this method con- 
verges to the correct two-body solution for small enough 
At. The rest of the paper will basically deal with the 
question of how small At needs to be in order for the 
method to produce meaningful results. 

To assess the convergence of the algorithm, we have nu- 
merically inverted the Laplace transform of {f), Eqs. Hll() 
and lO, for finite At using a publicly available Gaver- 
Wynn Rho algorithm ji^. In the following, we set units 
such that D = 1 and a = 1. The exact solutions are 
compared with averages obtained from simulation runs 
in two dimensions in Fig. ^ Each data point in the fig- 
ure corresponds to an average over 50000 runs starting 
from a given (xq, yo), each performing one step of size At 
as indicated in the different panels of the figure. 

De Michele's algorithm (circles) and Strating's algo- 
rithm (diamonds) give almost identical results for small 
enough At; only at unrealistically large At does one 
observe a difference due to 'missed' collisions in Strat- 
ing's algorithm. At this point, however, both algo- 



rithms already deviate significantly from the exact solu- 
tion (shown as a solid line) . The convergence to the right 
solution is good enough to give reasonable results already 
for At = 0.1, about one order of magnitude bi gger than 
what has been used in previous studies ITol lllj . For 
even smaller At, also the exact solution becomes indis- 
tinguishable from the one-dimensional (Ar)iD, the time- 
domain version of Eq. H1U|) : 

(^,.„.eta,^.--'.-.>_,„„,(_^) 

(13) 

(where po = ?'o — cr). 

Let us point out that this convergence is indeed closely 
related to the clastic-collision rule. In fact, from a 
Metropolis-MC scheme with Gaussian displacements , 
one gets the results shown as triangles in the figure, ap- 
proaching for small At its one-dimensional limit 



(Ar)iD.MC = e(po)^^e-''o/(4^'^*) . (14) 

Here we have dropped a term due to the possible 'tun- 
neling' of the particle across the excluded- volume region, 
which is exponentially small for At ^ cr^/ (4Z?). The two 
expressions, Eqs. (|13() and p4|l . differ in their leading- 
order terms by a factor of two. This highlights that the 
approach of the MC dynamics to true Brownian dynam- 
ics even as At ^ is less straightforward than one might 
suppose: the MC method neglects a term of ©(VAt) 
that is part of the finitc-time-step solution of the Smolu- 
chowski equation, whereas De Michele's algorithm has a 
leading error term of ©(At). A similar ar gum ent has 
been brought forward by Heyes and Braiika [43 for the 
case of smooth forces. There, however, one needed to 
look at the mean-squared displacement, (Ar^) to find 
differences at finite At, whereas here the disagreement 
sets in one level earlier. It also explains why the use of 
Monte Carlo to simulate BD needed an elaborate density- 
dependent extrapolation procedure to At = 0, using sev- 
eral simulation runs at different At > 

The Tao et al. ^3 algorithm is more difficult to assess. 
It replaces the deterministic reflection of the trajectory 
by a stochastic collision law and can for hard spheres be 
formulated as follows: assigning a fictive velocity to the 
particle from the random displacement Ar, one calculates 
the time < tc < At for which x{tc) = cr is reached, and 
then assigns a final position from new random displace- 
ments in accord with the remaining time At — tc- The 
procedure yields in one dimension 

PTao.iDl?", to + At|ro,to) = Go(r,to + At|ro,to) 

r-(a-ra)/At ^^^~v^ / {4DAt) 2g-^(r-o-) V(4DAi(l-r))) 

+ / \ = — dv, 

y_oo \/47rL»At ^47ri:>At(l - ry) 

(15) 

where rj = (cr— ro)/(vAt). The first term in the integral is 
the probability of assigning a fictive velocity v to the par- 
ticle such that a collision takes place at tc — {(J ~ ro)/v, 
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while the second term is the probabiHty to arrive at the 
final point r by a random displacement starting from 
r{tc) = cr and constrained to be directed away from the 
wall. Indeed, as rg a, we have rj —^ 0, and the sec- 
ond term in the above integral gives the correct image- 
point contribution, while the factor 2 cancels with the 
w-intcgration over the first term. Apart from this asymp- 
totic case, it is however difficult to evaluate pxao analyti- 
cally. Numerical evaluation shows it to be less accurate if 
(tq — a)/ At is not small but of intermediate magnitude. 
There is an ambiguity in the translation of the Tao-BD 
recipe described in Ref. when the coUiding objects 
have curved surfaces, regarding how the random post- 
collision velocities should be generated. Tao et al. state 
that they should be such that the particles separate ini- 
tially. If one implements this algorithm for hard spheres, 
the results for (Ar) arc almost identical to Dc Michele's 
ones for large At, as the plus symbols in Fig. show. 
For small At, a small but discernible difference always 
remains, owing to the randomness in the velocity reflec- 
tion. One can introduce a slight modification to the Tao- 
BD algorithm, in which the post-collision velocities arc 
only required to be such that the two colliding particles 
do not overlap at the end of the time step. If one does 
so, one can improve on the result for (A?^, as the crosses 
in Fig. n show: for small At, this algorithm behaves just 
like the original Tao-BD one, as expected, and for larger 
At, its results remain closer to the true solution. We did, 
however, not test this modification for the many-particle 
case: relaxing the criterion for the post-collision veloci- 
ties can in principle lead to the same secondary-overlap 
problems the Strating algorithm suffers from. 

Having discussed the one-step behavior of the algo- 
rithm, let us now look at the results after many steps 
M, i.e., at a time T = MAt large compared to At, 
M ^ 1. Here we perform a numerical comparison of 
the one-particle PDF in the two-dimensional case in the 
presence of fixed other spheres (corresponding to the rel- 
ative part of the PDF in the two-body problem). We re- 
strict the discussion to De Michele's algorithm now. As 
initial condition, let us fix a relative distance rg = I.Ict. 
Since we are interested in the PDF as a function oft, and 
the exact solution of this two-body case in Refs. |43. l4l| 
is given as a infinite scries in Laplace-transformed fre- 
quency s, we found it easier to solve for the PDF nu- 
merically. Using Crank's method [s^l on a 2-dimensional 
K X K periodic square mesh {K = 800) of length L 8, 
conditions of no flux are imposed on the circle; with 
such conditions Crank's method conserves the probabil- 
ity. Indicating with Sx = L/K, Crank's method is sta- 
ble for time steps D 5t <^ Sx^ ; the integration time step 
was chosen as 5t = lQ~'^5x'^ / D, checking that smaller 
time steps do not improve the accuracy of the numer- 
ical solution. The simulation data was obtained from 
averages over 200000 runs, then using data binning with 
dx = dy = 0.32 for both the simulation data and the 
numerically solved PDF. 

A visual inspection already yields some insight on the 
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FIG. 2: Contour plot comparing the probability distribu- 
tion (PDF) for the relative diffusion of two hard spheres (di- 
ameter a = 1, relative diffusion coefficient D = 1), at time 
T = 1/4. Circles (squares, diamonds) connected with dotted 
lines are the positions where the binned simulation-data PDF 
has reached its half-maximum (1/5, 1/50); filled symbols cor- 
respond to At = 0.01 (5000 steps), open symbols to At = 0.1 
(50 steps). Solid lines are the corresponding exact-PDF re- 
sults. A bin size of dx = dy = 0.32 was used. The grey circle 
indicates the position of the second particle (no-ffux boundary 
condition along the dashed line). 



quality of the agreement. To this end, we plot equal- 
probability contour lines corresponding to the half-, 1/5-, 
and 1 /50-maximum value of the PDF. The comparison of 
the exact result with De Michele's algorithm is shown in 
Fig. 121 There, solid lines indicate the exact results, while 
the different symbols represent the results of De Michele's 
algorithm with different time steps. At = 0.01 (corre- 
sponding to M = 5000 steps) and At = 0.1 (M = 50). 
The fluctuations visible in the figure are all well within 
the fluctuations expected for the number of runs used to 
obtain the average. 

While the above discussion refers to a two-dimensional 
test, the same features hold in three dimensions. There, a 
simple expression for the exact two-particle PDF can be 
obtained, if one starts from an angular-averaged initial 
distribution P(r,0|ro) = S{r — ro)/(47rrQ), viz. [sif : 



47rrrop(r,t|ro) 
1 



V47rDt 
1 



■ exp 



cxp 
Dt 



-{r - r^Y 



ADt 



+ exp 



erfc 




(16) 



(where Tq = 2(7 — tq). 

The reduced PDF, Eq. H16|) . is compared in Fig.|2|with 
the results from De Michele's algorithm. We show a cut 
using ro = l.Olcr and a fixed time step At = 0.1 in the 
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FIG. 3: Two-particle reduced PDF in three dimensions (with 
angular- averaged initial conditions), p{r,t\ro), Eq. 11611 . at 
times T — 0.1, 0.5, and 2.5, as a function of the scaled ini- 
tial separation (ro — a)/y/Dt, where ro = 1.01a (solid lines). 
Results from the De Michele's algorithm with At = 0.1 are 
shown as symbols; they correspond to M = 1, 5, and 25 sim- 
ulation time steps. 

algorithm. The final time T is varied through the num- 
ber of simulation time steps M. Note that this At, as 
judged from Fig.^ is already quite large, so that for only 
one simulation time step, a discrepancy is clearly visible. 
This however vanishes quite quickly; already for A/ w 5 
the deviations are minute, and for M w 25, they arc no 
longer identified. Generally, this indicates that the tran- 
sition from two to three dimension in the analysis of the 
algorithm poses no surprises. Furthermore, even if in the 
one-step comparison some differences remain for moder- 
ately large At, the results obtained after M 3> 1 steps 
can still be correct. 



IV. MANY PARTICLES CASE 

Wc now turn to a brief check of the many-body behav- 
ior of the algorithm. Here, analytical results are avail- 
able for the long-time diffusion coefficient in the low- 
density hmit [43 ■ One gets Dl = D{1 — 2(j)), where 
(p = (7r/6)ntT^ is the packing fraction, and n the number 
density of the hard-sphere system. Since this result is de- 
rived from the two-body PDF discussed above, it poses 
a genuine test of BD including hard-core exclusion. 

For the test of Dc Michele's algorithm, we used N = 
1000 particles at two different volume fractions ~ 0.05 
and (f) = 0.30. For each volume fraction several simula- 
tions at different time steps At were performed. Results 
for the mean-squared displacement (MSD) are shown as 
symbols in Fig. 01 normalized to the free diffusion asymp- 
tote, d{t) = (Ar^)/{6Dt). The long-time asymptote, 
d{t) = 1 — 20 is reached for Dt ^ 1, and this number does 
not depend on the time step At. The effect of the finite 
step size is clearly seen at short times. At ^ 0.05, where 
the ballistic sub-intervals in the algorithm lead to a MSD 
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FIG. 4; Mean-squared displacement (Ar^) of hard-sphere 
systems at packing fractions (j> — 0.05 (upper panel) and 
(j) — 0.30 (lower panel), normalized by the free diffusion 
limit, 6Dt. Symbols are simulations with De Michele's al- 
gorithm using 1000 particles with various step sizes At as 
indicated in the figure. The dashed horizontal line indi- 
cates the short-time asymptote for true Brownian dynamics, 
d{t) = (Ar^)/(6-Dt) — 1. The dot-dashed line corresponds 
to the long-time asymptote as evaluated in first order in the 
packing fraction, d{t) — 1 ~ 24>. 



quadratic in time, i.e., d{t) ^ t. In agreement with the 
discussion above (cf. Fig.O, the algorithm needs about 
C(10) steps in order to reproduce correctly the Brownian 
dynamics. One hence needs to reduce the time step to 
At = 0(10"'^) in order to recover a window in which the 
proper Brownian short-time dynamics is visible. Other- 
wise, the artificial ballistic d{t) crosses over directly to 
the correct long-time behavior, at least for realistic (not 
too big) choices of At. It is reassuring that the improper 
treatment of the short-time dynamics does not influence 
the convergence to the correct long-time dynamics, and 
does not even induce an effective time scale (or an effec- 
tive free diffusion coefficient). In other words, if one is 
not interested in very early times, a reasonably large At 
can be chosen to obtain the dynamics at t ^ At. This is 
a clear advantage over other schemes, in which an elab- 
orate At ^ extrapolation needed to be applied. For 
example, in the MC algorithm we find reasonable agree- 
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mcnt agreement for d(t) when Gaussian displacements 
of variance 5r^ = ©(O.Olcr) per Cartesian coordinate are 
chosen. Since in each MC step, a trial move for a sin- 
gle particle is made, this displacement size is connected 
to a time step by Sr^ = 2ND At 3. For our system, 
N = 1000, corresponding hence to a time step At ^ 10^ 
in our units. 

As the lower panel of Fig. 0] shows, the same conclu- 
sions can be drawn from a less dilute case, (j) — 0.3. Here, 
the first-order result quoted for Dl is no longer accurate. 
Still the Dl read off from the results of De Michele's al- 
gorithm is stable for At <C 0.1. At At = 0.1, a small 
deviation remains visible, and only at about At w 0.3 
(cf. the largest step size shown in the figure), the data 
show a major deviation also at long times. This data 
still resembles the value expected from the dilute case, 
d{t) = 1 — 2(f), which neglects three-body terms. The 
large step size thus violates the condition that the inter- 
vals At must deal with two-body 'events' only, demon- 
strating that the restriction At <^ dl^g/{4:D) for the max- 
imum tolerable stepsize indeed becomes the crucial one 
at moderate and high densities. 



V. CONCLUSIONS 

We have discussed schemes to integrate the Brownian 
motion of a many-body hard-sphere system, where the 
singular potential prevents the use of standard Brownian 
Dynamics techniques dealing with smooth forces. Em- 
phasis was placed on testing their correctness for small 
but finite integration timcsteps At. We have employed a 
set of tests based on the exactly known two-parti cle p rob- 
ability distribution functions for hard spheres ji^ . |4J|. 
These tests are sensitive to both the proper implemen- 
tation of free diffusion, i.e. the random force in the 
Langevin equation, Eq. , and to the treatment of hard- 
core 'collisions', i.e. the singular hard-sphere force in 
this equation. Hence they test both crucial ingredients 
to the HS-BD problem. 

It was shown that De Michele's ED-BD algorithm in- 
deed converges to the correct solution of Eq. ©. The 
algorithm works with a finite time step and performs a 
Newtonian dynamics simulation within each interval At, 
where the masses of the particles play the role of the 
inverse diffusion coefficients. Every At, random uncorre- 
lated fictive velocities are assigned to each particle, used 
as a book-keeping device to implement MC-like random 
moves without overlaps. This realizes the overdamped 
limit to the Langevin equation. Unlike a number of ear- 
lier schemes (cf. Ref. and citations therein), this ED- 
BD method avoids unphysical hard-sphere overlaps in 
any case. Methods based on soft-sphere approximations 
to the hard-sphere potential need to use very small inte- 
gration time steps, forced by the increasing steepness of 
the potential and the condition that potential forces vary 
little during any single time step. Monte Carlo methods, 
implementing hard-core exclusion but not flux reflection. 



reproduce Brownian dynamics strictly only for At 0. 
In contrast, by implementing the no-fiux boundary con- 
ditions at least approximately, De Michele's algorithm 
works reliably with significantly larger time steps. Their 
magnitude is not bound by the steepness of the potential, 
but only by the size of the hard spheres and their typical 
distance. The dynamical features of Eq. (jSJ are correctly 
reproduced after a small number of such time steps, and 
the method is stable in the sense that no drift or ef- 
fective rescaling is introduced to the long-time behavior 
when increasing the step size within the limits named 
above. This feature of De Michele's algorithm has been 
overlooked so far and is an improvement over methods 
requiring careful extrapolation to At ^ 0. 

As a simple extension of well-tested event-driven meth- 
ods for hard spheres, De Michele's algorithm shares their 
effectiveness (but also the problem of being difficult to 
parallelize). The same holds for the method recently de- 
veloped by Tao et al. although the latter is some- 
what more expensive in terms of computing time: there, 
considerably more random numbers need to be drawn 
(at least one more per collision and dimension), while 
Dc Michele's algorithm uses simpler velocity-reflection 
laws. We found no considerable difference concerning 
the results between the two methods. 

De Michele's algorithm can easily be extended from the 
equilibrium case tested here. We have discussed in par- 
ticular how to include linear shear flows (of in principle 
arbitrary magnitude). At large shear rate, one needs to 
take into account two modifications, namely a determin- 
istic drift and a distortion of the fictive- velocity distribu- 
tion. The latter has so far been ignored in the discussion 
of BD-HS algorithms. Further extensions such as the 
inclusion of finite inertial terms, or of finite stepwise in- 
teractions, will be discussed in subsequent publications. 
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APPENDIX A: AVERAGE RELATIVE 
DISPLACEMENTS 

We calculate the average relative displacement of two 
Brownian hard spheres of diameter a; fixing the rela- 
tive diffusion coefficient to be D. This quantity is most 
conveniently calculated employing Laplace-transformed 
time (frequency s = Dq^) and partial integration of the 
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diffusion equation, viz. in ID 

/•OO 

(r(s))iD = {D/q^) / rdrp{r, s\ro) dx 



iD/q^)p{a,s\ro)^{l/q')e- 



-q{ra-u) 



(Al) 



27r 



d^[(Tp(cr, 0, s|ro,0o)cos(/)] . (A2) 



Replacing cos0 by sin^ yields {y{s)). Inserting Eq. lO, 
p{r,(f>,s\ro,(j>o) = (l/D)X;™Cxp[im((/) - 0o)]Pm('',s|ro), 
one notes that only the m = ±1 terms contribute to 
(x) and (y). Furthermore, for r = a, the square brackets 
simplify by virtue of the Wronskian of the modified Bessel 
fimctions, K'„^{z)Im{z) — Km{z)I'm{z) ~ 1/z, so that 



where p{x,s\xq) is the two-particle PDF in the Laplace 
domain (see Ref. j44j ) and boundary terms vanish due to 
the no-flux boundary condition. 

A similar result holds in 2D: Note that for the Carte- 
sian coordinate x, 

/•C30 /•27r 

{q^lD){x{s))2^ - / dr d0 [rdr{rdrp) + S^p] cos(A 



Pi(ct, s|ro) = 



1 K^{qro) 
q<T K[{qa) 



(A3) 



Summarizing these expressions, one arrives at Eq. Hll|l . 
The 3D result, Eq. H12|l is obtained along the same lines. 
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